#include "galois.h"
#include <stdlib.h>
#include <stdio.h>
#include <vector>
#include <iostream>

#define rgw32_mask(v) ((v) & 0x80000000)

/* type: -1 == choose by w, 0 == table, 1 == log, 2 == shift, 3 == split (32 only) */
/* Hard limit on table? */
Galois *galois_init(int w, int type)
{
	Galois *gc;
	if (w < 0 || w > 32) {
		fprintf(stderr, "No implementation for w = %d\n", w);
		exit(1);
	}
	if (type == -1) {
		if (w > 0 && w < 10) {
			gc = new Galois_Table(w);
			return gc;
		} else if (w < 23) {
			gc = new Galois_Log(w);
			return gc;
		} else if (w < 32) {
			gc = new Galois_Shift(w);
			return gc;
		} else if (w == 32) {
			gc = new Galois_Split(w);
			return gc;
		}

	} else if (type == 0) {
		/*if (w > 9) {
			fprintf(stderr, "W value too large for Tables, w = %d\n", w);
			exit(1);
		} */ /* Limit? */;
		gc = new Galois_Table(w);
		return gc;

	} else if (type == 1) {
		gc = new Galois_Log(w);
		return gc;

	} else if (type == 2) {
		gc = new Galois_Shift(w);
		return gc;

	} else if (type == 3) {
		if (w != 32) {
			cerr << "Split table construction only for w = 32\n";
			exit(1);
		}
		gc = new Galois_Split(w);
		return gc;
	}
}

/***** TABLE METHOD DEFINITIONS *****/
Galois_Table::Galois_Table(int wtmp) 
{
	w = wtmp;
	if(w <= 8) {
		if (galois_create_mult_tables(w, &galois_w08_mult_table, &galois_w08_div_table) < 0) {
			fprintf(stderr, "Error: create_galois_mult_tables\n");
			exit(1);
		}
	} else if (w <= 16) {
		if (galois_create_mult_tables(w, &galois_w16_mult_table, &galois_w16_div_table) < 0) {
			fprintf(stderr, "Error: create_galois_mult_tables\n");
			exit(1);
		}
	} else {
		fprintf(stderr, "No table implementation for that value of w\n");
		exit(1);
	}
}

int Galois_Table::galois_single_multiply(int a, int b)
{
	if (w < 9) return galois_w08_mult_table[(a << w) | b];
	else return galois_w16_mult_table[(a << w) | b];
}

int Galois_Table::galois_single_divide(int a, int b)
{
	if (w < 9) return galois_w08_div_table[(a << w) | b];
	else return galois_w16_div_table[(a << w) | b];
}

int Galois_Table::galois_inverse(int a)
{
	return galois_single_divide(1, a);
}

void Galois_Table::galois_w08_region_multiply(unsigned char *region, int multby, int nbytes,
		unsigned char *r2, int add)
{
	unsigned char *ur1, *ur2, *cp;
	unsigned char prod;
	int i, srow, j;
	uint64_t l, *lp2;
	unsigned char *lp;
	int sol;

	ur1 = (unsigned char *) region;
	ur2 = (r2 == NULL) ? ur1 : (unsigned char *) r2;

	if (w != 8) {
		fprintf(stderr, "Galois_Table::galois_w08_region_multiply called with w == %d\n", w);
		exit(1);
	}

	/* With galois w==08, we multiply single characters (int value 0 - 255).
	   fixed width uint64_t are used as the packet size, but it is guaranteed
	   to be that size at this point */
	srow = multby * (1 << 8);
	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			prod = galois_w08_mult_table[srow+ur1[i]];
			ur2[i] = prod;
		}
	} else {
		sol = sizeof(uint64_t);
		lp2 = &l;
		lp = (unsigned char *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = ur2+i;
			lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				prod = galois_w08_mult_table[srow+ur1[i+j]];
				lp[j] = prod;
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	return;

}

void Galois_Table::galois_w16_region_multiply(unsigned char *region, int multby, int nbytes,
		unsigned char *r2, int add)
{
	/*
	   uint16_t *ur1, *ur2, *cp;
	   uint16_t prod;
	   int i, srow, j;
	   uint64_t l, *lp2, *lptop;
	   uint16_t *lp;
	   int sol;

	   ur1 = (uint16_t *) region;
	   ur2 = (r2 == NULL) ? ur1 : (uint16_t *) r2;
	   nbytes /= sizeof(uint16_t);

	   if (multby == 0) {
	   if (!add) {
	   lp2 = (uint64_t *) ur2;
	   ur2 += nbytes;
	   lptop = (uint64_t *) ur2;
	   while (lp2 < lptop) { *lp2 = 0; lp2++; };
	   }
	   return;
	   }

	//With galois w==08, we multiply single characters (int value 0 - 255).
	//fixed width uint64_t are used as the packet size, but it is guaranteed
	//to be that size at this point
	//srow = multby * (1 << 16);
	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			prod = galois_w16_mult_tables[8][srow+ur1[i]];
			ur2[i] = prod;
		}
	} else {
		sol = sizeof(uint64_t)/2;
		lp2 = &l;
		lp = (uint16_t *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = ur2+i;
			lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				prod = galois_w16_mult_tables[8][srow+ur1[i+j]];
				lp[j] = prod;
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	return;
	*/
}

void Galois_Table::galois_w32_region_multiply(unsigned char *region, int multby, int nbytes,
		unsigned char *r2, int add) {}


/***** LOG METHOD DEFINITIONS *****/
Galois_Log::Galois_Log(int wtmp) {
	w = wtmp;
	if (w <= 8) {
		if (galois_create_log_tables(w, &galois_w08_log_table, &galois_w08_ilog_table) < 0) {
			fprintf(stderr, "Error: create_galois_log_tables\n");
			exit(1);
		}
	} else if (w<=16) {
		if (galois_create_log_tables(w, &galois_w16_log_table, &galois_w16_ilog_table) < 0) {
			fprintf(stderr, "Error: create_galois_log_tables\n");
			exit(1);
		}
	} else if (w <= 32) {
		if (galois_create_log_tables(w, &galois_w32_log_table, &galois_w32_ilog_table) < 0) {
			fprintf(stderr, "Error: create_galois_log_tables\n");
			exit(1);
		}
	} else {
		fprintf(stderr, "Error: no log implementation for that value of w\n");
		exit(1);
	}
}

int Galois_Log::galois_single_multiply(int a, int b)
{
	int sum_j, z;
	if (a == 0 || b == 0) return 0;
	if (w < 9) {
		sum_j = (int) galois_w08_log_table[a] + (int) galois_w08_log_table[b];
		z = (int) galois_w08_ilog_table[sum_j];
		return z;
	} else if (w < 17) {
		sum_j = (int) galois_w16_log_table[a] + (int) galois_w16_log_table[b];
		z = (int) galois_w16_ilog_table[sum_j];
		return z;
	} else if (w < 33) {
		sum_j = (int) galois_w32_log_table[a] + (int) galois_w32_log_table[b];
		z = (int) galois_w32_ilog_table[sum_j];
		return z;
	}
}

int Galois_Log::galois_single_divide(int a, int b)
{
	int sum_j, z;
	if (b == 0) return -1;
	if (a == 0) return 0;
	if (w < 9) {
		sum_j = (int) galois_w08_log_table[a] - (int) galois_w08_log_table[b];
		z = (int) galois_w08_ilog_table[sum_j];
		return z;
	} else if (w < 17) {
		sum_j = (int) galois_w16_log_table[a] - (int) galois_w16_log_table[b];
		z = (int) galois_w16_ilog_table[sum_j];
		return z;
	} else if (w < 33) {
		sum_j = (int) galois_w32_log_table[a] - (int) galois_w32_log_table[b];
		z = (int) galois_w32_ilog_table[sum_j];
		return z;
	}
}

int Galois_Log::galois_inverse(int a)
{
	return galois_single_divide(1, a);
}

void Galois_Log::galois_w08_region_multiply(unsigned char *region, int multby, int nbytes,
		unsigned char *r2, int add)
{
	unsigned char *ur1, *ur2, *cp;
	int prod;
	int i, log1, j, log2;
	uint64_t l, *lp2, *lptop;
	unsigned char *lp;
	int sol;

	ur1 = (unsigned char *) region;
	ur2 = (r2 == NULL) ? ur1 : (unsigned char *) r2;

	if (multby == 0) {
		if (!add) {
			lp2 = (uint64_t *) ur2;
			ur2 += nbytes;
			lptop = (uint64_t *) ur2;
			while (lp2 < lptop) { *lp2 = 0; lp2++; }
		}
		return;
	}

	log1 = (int) galois_w08_log_table[multby];

	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			if (ur1[i] == 0) {
				ur2[i] = 0;
			} else {
				prod = (int) galois_w08_log_table[ur1[i]] + log1;
				ur2[i] = galois_w08_ilog_table[prod];
			}
		}
	} else {
		sol = sizeof(uint64_t);
		lp2 = &l;
		lp = (unsigned char *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = ur2+i;
			lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				if (ur1[i+j] == 0) {
					lp[j] = 0;
				} else {
					log2 = (int) galois_w08_log_table[ur1[i+j]];
					prod = log2 + log1;
					lp[j] = galois_w08_ilog_table[prod];
				}
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	return;

}

void Galois_Log::galois_w16_region_multiply(unsigned char *region, int multby, int nbytes,
		unsigned char *r2, int add)
{
	uint16_t *ur1, *ur2, *cp;
	int  prod;
	int i, log1, j, log2;
	uint64_t l, *lp2, *lptop;
	uint16_t *lp;
	int sol;

	ur1 = (uint16_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint16_t *) r2;
	nbytes /= sizeof(uint16_t);

	if (multby == 0) {
		if (!add) {
			lp2 = (uint64_t *) ur2;
			ur2 += nbytes;
			lptop = (uint64_t *) ur2;
			while (lp2 < lptop) { *lp2 = 0; lp2++; }
		}
		return;
	}

	log1 = galois_w16_log_table[multby];

	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			if (ur1[i] == 0) {
				ur2[i] = 0;
			} else {
				prod = galois_w16_log_table[ur1[i]] + log1;
				ur2[i] = galois_w16_ilog_table[prod];
			}
		}
	} else {
		sol = sizeof(uint64_t)/2;
		lp2 = &l;
		lp = (uint16_t *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = ur2+i;
			lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				if (ur1[i+j] == 0) {
					lp[j] = 0;
				} else {
					log2 = galois_w16_log_table[ur1[i+j]];
					prod = log2 + log1;
					lp[j] = galois_w16_ilog_table[prod];
				}
			}
			*lp2 = (*lp2) ^ l;
		}
	}
	return;

}

void Galois_Log::galois_w32_region_multiply(unsigned char *region, int multby, int nbytes,
		unsigned char *r2, int add)
{
}

/***** SHIFT METHOD DEFINITIONS ******/
Galois_Shift::Galois_Shift(int wtmp)
{
	w = wtmp;
}

int Galois_Shift::galois_single_multiply(int x, int y)
{
	int prod;
	int i, j, ind;
	int k;
	int scratch[33];

	int prim_poly[33] =
	{ 0,
		/*  1 */     1,
		/*  2 */    07,
		/*  3 */    013,
		/*  4 */    023,
		/*  5 */    045,
		/*  6 */    0103,
		/*  7 */    0211,
		/*  8 */    0435,
		/*  9 */    01021,
		/* 10 */    02011,
		/* 11 */    04005,
		/* 12 */    010123,
		/* 13 */    020033,
		/* 14 */    042103,
		/* 15 */    0100003,
		/* 16 */    0210013,
		/* 17 */    0400011,
		/* 18 */    01000201,
		/* 19 */    02000047,
		/* 20 */    04000011,
		/* 21 */    010000005,
		/* 22 */    020000003,
		/* 23 */    040000041,
		/* 24 */    0100000207,
		/* 25 */    0200000011,
		/* 26 */    0400000107,
		/* 27 */    01000000047,
		/* 28 */    02000000011,
		/* 29 */    04000000005,
		/* 30 */    010040000007,
		/* 31 */    020000000011,
		/* 32 */    00020000007 };  /* Really 40020000007, but we're omitting the high order bit */

	int nwm1[33] = { 0, (1 << 1)-1, (1 << 2)-1, (1 << 3)-1, (1<< 4)-1,
		(1 << 5)-1, (1 << 6)-1, (1 << 7)-1, (1 << 8)-1, (1 << 9)-1, (1 << 10)-1,
		(1 << 11)-1, (1 << 12)-1, (1 << 13)-1, (1 << 14)-1, (1 << 15)-1, (1 << 16)-1,
		(1 << 17)-1, (1 << 18)-1, (1 << 19)-1, (1 << 20)-1, (1 << 21)-1, (1 << 22)-1,
		(1 << 23)-1, (1 << 24)-1, (1 << 25)-1, (1 << 26)-1, (1 << 27)-1, (1 << 28)-1,
		(1 << 29)-1, (1 << 30)-1, 0x7fffffff, 0xffffffff };

	/* set values of scratch to be log values, starting at y */
	prod = 0;
	for (i = 0; i < w; i++) {
		scratch[i] = y;
		if (y & (1 << (w-1))) {
			y = y << 1;
			y = (y ^ prim_poly[w]) & nwm1[w];
		} else {
			y = y << 1;
		}
	}
	/* do magic */
	for (i = 0; i < w; i++) {
		ind = (1 << i);
		if (ind & x) {
			j = 1;
			for (k = 0; k < w; k++) {
				prod = prod ^ (j & scratch[i]);
				j = (j << 1);
			}
		}
	}
	return prod;
}

int Galois_Shift::galois_single_divide(int a, int b)
{
	int inverse;

	if (b == 0) return -1;
	if (a == 0) return 0;
	inverse = galois_inverse(b);
	return galois_single_multiply(a, inverse);
}

int Galois_Shift::galois_inverse(int y)
{
	int mat[1024], mat2[32];
	int inv[1024], inv2[32];
	int ind, i, j, k, prod;

	int prim_poly[33] =
	{ 0,
		/*  1 */     1,
		/*  2 */    07,
		/*  3 */    013,
		/*  4 */    023,
		/*  5 */    045,
		/*  6 */    0103,
		/*  7 */    0211,
		/*  8 */    0435,
		/*  9 */    01021,
		/* 10 */    02011,
		/* 11 */    04005,
		/* 12 */    010123,
		/* 13 */    020033,
		/* 14 */    042103,
		/* 15 */    0100003,
		/* 16 */    0210013,
		/* 17 */    0400011,
		/* 18 */    01000201,
		/* 19 */    02000047,
		/* 20 */    04000011,
		/* 21 */    010000005,
		/* 22 */    020000003,
		/* 23 */    040000041,
		/* 24 */    0100000207,
		/* 25 */    0200000011,
		/* 26 */    0400000107,
		/* 27 */    01000000047,
		/* 28 */    02000000011,
		/* 29 */    04000000005,
		/* 30 */    010040000007,
		/* 31 */    020000000011,
		/* 32 */    00020000007 };  /* Really 40020000007, but we're omitting the high order bit */

	int nw[33] = { 0, (1 << 1), (1 << 2), (1 << 3), (1 << 4),
		(1 << 5), (1 << 6), (1 << 7), (1 << 8), (1 << 9), (1 << 10),
		(1 << 11), (1 << 12), (1 << 13), (1 << 14), (1 << 15), (1 << 16),
		(1 << 17), (1 << 18), (1 << 19), (1 << 20), (1 << 21), (1 << 22),
		(1 << 23), (1 << 24), (1 << 25), (1 << 26), (1 << 27), (1 << 28),
		(1 << 29), (1 << 30), (1 << 31), -1 };


	int nwm1[33] = { 0, (1 << 1)-1, (1 << 2)-1, (1 << 3)-1, (1<< 4)-1,
		(1 << 5)-1, (1 << 6)-1, (1 << 7)-1, (1 << 8)-1, (1 << 9)-1, (1 << 10)-1,
		(1 << 11)-1, (1 << 12)-1, (1 << 13)-1, (1 << 14)-1, (1 << 15)-1, (1 << 16)-1,
		(1 << 17)-1, (1 << 18)-1, (1 << 19)-1, (1 << 20)-1, (1 << 21)-1, (1 << 22)-1,
		(1 << 23)-1, (1 << 24)-1, (1 << 25)-1, (1 << 26)-1, (1 << 27)-1, (1 << 28)-1,
		(1 << 29)-1, (1 << 30)-1, 0x7fffffff, 0xffffffff };

	for (i = 0; i < w; i++) {
		mat2[i] = y;

		if (y & nw[w-1]) {
			y = y << 1;
			y = (y ^ prim_poly[w]) & nwm1[w];
		} else {
			y = y << 1;
		}
	}

	galois_invert_binary_matrix(mat2, inv2, w);

	return inv2[0];
}

void Galois_Shift::galois_w08_region_multiply(unsigned char *region, int multby, int nbytes,
		unsigned char *r2, int add) { }
void Galois_Shift::galois_w16_region_multiply(unsigned char *region, int multby, int nbytes,
		unsigned char *r2, int add) { }
void Galois_Shift::galois_w32_region_multiply(unsigned char *region, int multby, int nbytes,
		unsigned char *r2, int add)
{
	uint16_t *ur1, *ur2, *cp;
	uint16_t prod;
	int i, srow, j;
	uint64_t l, *lp2, *lptop;
	uint16_t *lp;
	int sol;

	ur1 = (uint16_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint16_t *) r2;
	nbytes /= sizeof(uint16_t);

	if (multby == 0) {
		if (!add) {
			lp2 = (uint64_t *) ur2;
			ur2 += nbytes;
			lptop = (uint64_t *) ur2;
			while (lp2 < lptop) { *lp2 = 0; lp2++; };
		}
		return;
	}

	/* With galois w==08, we multiply single characters (int value 0 - 255).
	   fixed width uint64_t are used as the packet size, but it is guaranteed
	   to be that size at this point */
	if (r2 == NULL || !add) {
		for (i = 0; i < nbytes; i++) {
			ur2[i] = galois_single_multiply(multby, ur1[i]);
		}
	} else {
		sol = sizeof(uint64_t)/2;
		lp2 = &l;
		lp = (uint16_t *) lp2;
		for (i = 0; i < nbytes; i += sol) {
			cp = ur2+i;                        lp2 = (uint64_t *) cp;
			for (j = 0; j < sol; j++) {
				lp[j] = galois_single_multiply(multby, ur1[i+j]);
			}
			*lp2 = (*lp2) ^ l;
		}
	}        return;

}

/***** SPLIT METHOD DEFINITIONS ******/
Galois_Split::Galois_Split(int wtmp)
{
	w = wtmp;
	galois_split_w08_tables = new uint32_t * [7];
	if (wtmp != 32) {
		fprintf(stderr, "Split tables only implemented for w = 32\n");
		exit(1);
	}

	if (galois_create_split_w08_tables(galois_split_w08_tables) < 0) {
		fprintf(stderr, "Error: create_galois_split_w08_tables\n");
		exit(1);
	}
}

int Galois_Split::galois_single_multiply(int x, int y)
{
	int i, j, a, b, accumulator, i8, j8;

	accumulator = 0;

	i8 = 0;
	for (i = 0; i < 4; i++) {
		a = (((x >> i8) & 255) << 8);
		j8 = 0;
		for (j = 0; j < 4; j++) {
			b = ((y >> j8) & 255);
			accumulator ^= galois_split_w08_tables[i+j][a|b];
			j8 += 8;
		}
		i8 += 8;
	}
	return accumulator;
}

int Galois_Split::galois_single_divide(int a, int b)
{
	int x;

	if (b == 0) return -1;
	if (a == 0) return 0;

	x = galois_inverse(b);
	return galois_single_multiply(a, x);
}

int Galois_Split::galois_inverse(int y)
{
	int mat[1024], mat2[32];
	int inv[1024], inv2[32];
	int ind, i, j, k, prod;

	int prim_poly[33] =
	{ 0,
		/*  1 */     1,
		/*  2 */    07,
		/*  3 */    013,
		/*  4 */    023,
		/*  5 */    045,
		/*  6 */    0103,
		/*  7 */    0211,
		/*  8 */    0435,
		/*  9 */    01021,
		/* 10 */    02011,
		/* 11 */    04005,
		/* 12 */    010123,
		/* 13 */    020033,
		/* 14 */    042103,
		/* 15 */    0100003,
		/* 16 */    0210013,
		/* 17 */    0400011,
		/* 18 */    01000201,
		/* 19 */    02000047,
		/* 20 */    04000011,
		/* 21 */    010000005,
		/* 22 */    020000003,
		/* 23 */    040000041,
		/* 24 */    0100000207,
		/* 25 */    0200000011,
		/* 26 */    0400000107,
		/* 27 */    01000000047,
		/* 28 */    02000000011,
		/* 29 */    04000000005,
		/* 30 */    010040000007,
		/* 31 */    020000000011,
		/* 32 */    00020000007 };  /* Really 40020000007, but we're omitting the high order bit */

	int nw[33] = { 0, (1 << 1), (1 << 2), (1 << 3), (1 << 4),
		(1 << 5), (1 << 6), (1 << 7), (1 << 8), (1 << 9), (1 << 10),
		(1 << 11), (1 << 12), (1 << 13), (1 << 14), (1 << 15), (1 << 16),
		(1 << 17), (1 << 18), (1 << 19), (1 << 20), (1 << 21), (1 << 22),
		(1 << 23), (1 << 24), (1 << 25), (1 << 26), (1 << 27), (1 << 28),
		(1 << 29), (1 << 30), (1 << 31), -1 };


	int nwm1[33] = { 0, (1 << 1)-1, (1 << 2)-1, (1 << 3)-1, (1<< 4)-1,
		(1 << 5)-1, (1 << 6)-1, (1 << 7)-1, (1 << 8)-1, (1 << 9)-1, (1 << 10)-1,
		(1 << 11)-1, (1 << 12)-1, (1 << 13)-1, (1 << 14)-1, (1 << 15)-1, (1 << 16)-1,
		(1 << 17)-1, (1 << 18)-1, (1 << 19)-1, (1 << 20)-1, (1 << 21)-1, (1 << 22)-1,
		(1 << 23)-1, (1 << 24)-1, (1 << 25)-1, (1 << 26)-1, (1 << 27)-1, (1 << 28)-1,
		(1 << 29)-1, (1 << 30)-1, 0x7fffffff, 0xffffffff };


	for (i = 0; i < w; i++) {
		mat2[i] = y;

		if (y & nw[w-1]) {
			y = y << 1;
			y = (y ^ prim_poly[w]) & (nwm1[w]);
		} else {
			y = y << 1;
		}
	}

	galois_invert_binary_matrix(mat2, inv2, w);

	return inv2[0];
}

void Galois_Split::galois_w08_region_multiply(unsigned char *region, int multby, int nbytes,
		unsigned char *r2, int add) {}
void Galois_Split::galois_w16_region_multiply(unsigned char *region, int multby, int nbytes,
		unsigned char *r2, int add) {}
void Galois_Split::galois_w32_region_multiply(unsigned char *region, int multby, int nbytes,
		unsigned char *r2, int add)
{
	uint32_t *ur1, *ur2, *cp, *ur2top;
	uint64_t *lp2, *lptop;
	int i, j, a, b, accumulator, i8, j8, k;
	int acache[4];

	ur1 = (uint32_t *) region;
	ur2 = (r2 == NULL) ? ur1 : (uint32_t *) r2;
	nbytes /= sizeof(uint32_t);
	ur2top = ur2 + nbytes;

	/* If we're overwriting r2, then we can't do better than just calling split_multiply.
	   We'll inline it here to save on the procedure call overhead */

	i8 = 0;
	for (i = 0; i < 4; i++) {
		acache[i] = (((multby >> i8) & 255) << 8);
		i8 += 8;
	}

	/* This works */
	if (!add) {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			for (i = 0; i < 4; i++) {
				a = acache[i];
				j8 = 0;
				for (j = 0; j < 4; j++) {
					b = ((ur1[k] >> j8) & 255);
					accumulator ^= galois_split_w08_tables[i+j][a|b];
					j8 += 8;
				}
			}
			ur2[k] = accumulator;
		}
	} else {
		for (k = 0; k < nbytes; k++) {
			accumulator = 0;
			for (i = 0; i < 4; i++) {
				a = acache[i];
				j8 = 0;
				for (j = 0; j < 4; j++) {
					b = ((ur1[k] >> j8) & 255);
					accumulator ^= galois_split_w08_tables[i+j][a|b];
					j8 += 8;
				}
			}
			ur2[k] = (ur2[k] ^ accumulator);
		}
	}
	return;
}

int galois_create_mult_tables(int w, unsigned char **galois_mult_table, unsigned char **galois_div_table)
{
	int j, x, y, logx;
	Galois_Log *gc;
	unsigned int prim_poly[9] =
	{ 0,
		/*  1 */     1,
		/*  2 */    07,
		/*  3 */    013,
		/*  4 */    023,
		/*  5 */    045,
		/*  6 */    0103,
		/*  7 */    0211,
		/*  8 */    0435,
	};

	unsigned int nw[9] = { 0, (1 << 1), (1 << 2), (1 << 3), (1 << 4),
		(1 << 5), (1 << 6), (1 << 7), (1 << 8) };

	unsigned int nwm1[9] = { 0, (1 << 1)-1, (1 << 2)-1, (1 << 3)-1, (1<< 4)-1,
		(1 << 5)-1, (1 << 6)-1, (1 << 7)-1, (1 << 8)-1 };

	try {
		galois_mult_table[0] = new unsigned char [nw[w] * nw[w]];
		galois_div_table[0] = new unsigned char [nw[w] * nw[w]];
	} catch (bad_alloc &e) {
		if (galois_mult_table[0] != NULL) delete galois_mult_table[0];
		if (galois_div_table[0] != NULL) delete galois_div_table[0];
		galois_mult_table[0] = NULL;
		galois_div_table[0] = NULL;
		return -1;
	}

	gc = new Galois_Log(w);
	if (gc == NULL) {
		delete galois_mult_table[0];
		delete galois_div_table[0];
		return -1;
	}

	/* Set mult/div tables for x = 0 */
	j = 0;
	galois_mult_table[0][j] = (unsigned char) 0;   /* y = 0 */
	galois_div_table[0][j] = (unsigned char) -1;
	j++;
	for (y = 1; y < nw[w]; y++) {   /* y > 0 */
		galois_mult_table[0][j] = (unsigned char) 0;
		galois_div_table[0][j] = (unsigned char) 0;
		j++;
	}

	for (x = 1; x < nw[w]; x++) {  /* x > 0 */
		galois_mult_table[0][j] = (unsigned char) 0; /* y = 0 */
		galois_div_table[0][j] = (unsigned char) -1;
		j++;
		logx = gc->galois_w08_log_table[x];
		for (y = 1; y < nw[w]; y++) {  /* y > 0 */
			galois_mult_table[0][j] = gc->galois_w08_ilog_table[logx+gc->galois_w08_log_table[y]];
			galois_div_table[0][j] = gc->galois_w08_ilog_table[logx-gc->galois_w08_log_table[y]];
			j++;
		}
	}
	return 0;

}

int galois_create_mult_tables(int w, uint16_t **galois_mult_table, uint16_t **galois_div_table)
{
	int j, x, y, logx;
	Galois_Log *gc;
	int prim_poly[17] ={ 0,
		/*  1 */     1,
		/*  2 */    07,
		/*  3 */    013,
		/*  4 */    023,
		/*  5 */    045,
		/*  6 */    0103,
		/*  7 */    0211,
		/*  8 */    0435,
		/*  9 */    01021,
		/* 10 */    02011,
		/* 11 */    04005,
		/* 12 */    010123,
		/* 13 */    020033,
		/* 14 */    042103,
		/* 15 */    0100003,
		/* 16 */    0210013 };

	int nw[17] = { 0, (1 << 1), (1 << 2), (1 << 3), (1 << 4),
		(1 << 5), (1 << 6), (1 << 7), (1 << 8), (1 << 9), (1 << 10),
		(1 << 11), (1 << 12), (1 << 13), (1 << 14), (1 << 15), (1 << 16) };


	int nwm1[17] = { 0, (1 << 1)-1, (1 << 2)-1, (1 << 3)-1, (1<< 4)-1,
		(1 << 5)-1, (1 << 6)-1, (1 << 7)-1, (1 << 8)-1, (1 << 9)-1, (1 << 10)-1,
		(1 << 11)-1, (1 << 12)-1, (1 << 13)-1, (1 << 14)-1, (1 << 15)-1, (1 << 16)-1 };

	try {
		galois_mult_table[0] = new uint16_t [nw[w] * nw[w]];
		galois_div_table[0] = new uint16_t [nw[w] * nw[w]];
	} catch (bad_alloc &e) {
		if (galois_mult_table[0] != NULL) delete galois_mult_table[0];
		if (galois_div_table[0] != NULL) delete galois_div_table[0];
		galois_mult_table[0] = NULL;
		galois_div_table[0] = NULL;
		return -1;
	}

	gc = new Galois_Log(w);
	if (gc == NULL) {
		delete galois_mult_table[0];
		delete galois_div_table[0];
		return -1;
	}

	/* Set mult/div tables for x = 0 */
	j = 0;
	galois_mult_table[0][j] = (uint16_t) 0;   /* y = 0 */
	galois_div_table[0][j] = (uint16_t) -1;
	j++;
	for (y = 1; y < nw[w]; y++) {   /* y > 0 */
		galois_mult_table[0][j] = (uint16_t) 0;
		galois_div_table[0][j] = (uint16_t) 0;
		j++;
	}

	for (x = 1; x < nw[w]; x++) {  /* x > 0 */
		galois_mult_table[0][j] = (uint16_t) 0; /* y = 0 */
		galois_div_table[0][j] = (uint16_t) -1;
		j++;
		logx = gc->galois_w16_log_table[x];
		for (y = 1; y < nw[w]; y++) {  /* y > 0 */
			galois_mult_table[0][j] = gc->galois_w16_ilog_table[logx+gc->galois_w16_log_table[y]];
			galois_div_table[0][j] = gc->galois_w16_ilog_table[logx-gc->galois_w16_log_table[y]];
			j++;
		}
	}
	return 0;
}

int galois_create_log_tables(int w, unsigned char **log_table, unsigned char **ilog_table)
{
	int i, j, b;
	int prim_poly[9] =
	{ 0,
		/*  1 */     1,
		/*  2 */    07,
		/*  3 */    013,
		/*  4 */    023,
		/*  5 */    045,
		/*  6 */    0103,
		/*  7 */    0211,
		/*  8 */    0435,
	};

	int nw[9] = { 0, (1 << 1), (1 << 2), (1 << 3), (1 << 4),
		(1 << 5), (1 << 6), (1 << 7), (1 << 8) }; 

	int nwm1[9] = { 0, (1 << 1)-1, (1 << 2)-1, (1 << 3)-1, (1<< 4)-1,
		(1 << 5)-1, (1 << 6)-1, (1 << 7)-1, (1 << 8)-1 };

	try {
		if (log_table[0] == NULL) log_table[0] = new unsigned char [nw[w]];
		if (ilog_table[0] == NULL) ilog_table[0] = new unsigned char [nw[w]*3];
	} catch (bad_alloc &e) {
		cerr << e.what() << endl;
		if (log_table[0] != NULL) delete log_table[0];
		if (ilog_table[0] != NULL) delete ilog_table[0];
		log_table[0] = NULL;
		ilog_table[0] = NULL;
		return -1;
	}

	for (j  = 0; j < nw[w]; j++) {
		log_table[0][j] = (unsigned char) nwm1[w];
		ilog_table[0][j] =  (unsigned char) 0;
	}

	b = 1;
	for (j = 0; j < nwm1[w]; j++) {
		if (log_table[0][b] != (unsigned char) nwm1[w]) {
			fprintf(stderr, "Galois_create_log_tables Error: j=%d, b=%d, B->J[b]=%d, J->B[j]=%d (0%o)\n",
					j, b, log_table[0][b], ilog_table[0][j], (b << 1) ^ prim_poly[w]);
			exit(1);
		}
		log_table[0][b] = (unsigned char) j;
		ilog_table[0][j] = (unsigned char) b;
		b <<= 1;
		if (b & nw[w]) b = (b ^ prim_poly[w]) & nwm1[w];
	}

	for (j = 0; j < nwm1[w]; j++) {
		ilog_table[0][j+nwm1[w]] = ilog_table[0][j];
		ilog_table[0][j+nwm1[w]*2] = ilog_table[0][j];
	}
	ilog_table[0] += nwm1[w];
	return 0;
}

int galois_create_log_tables(int w, uint16_t **log_table, uint16_t **ilog_table)
{
	int i, j, b;

	int prim_poly[17] ={ 0,
		/*  1 */     1,
		/*  2 */    07,
		/*  3 */    013,
		/*  4 */    023,
		/*  5 */    045,
		/*  6 */    0103,
		/*  7 */    0211,
		/*  8 */    0435,
		/*  9 */    01021,
		/* 10 */    02011,
		/* 11 */    04005,
		/* 12 */    010123,
		/* 13 */    020033,
		/* 14 */    042103,
		/* 15 */    0100003,
		/* 16 */    0210013 };

	int nw[17] = { 0, (1 << 1), (1 << 2), (1 << 3), (1 << 4),
		(1 << 5), (1 << 6), (1 << 7), (1 << 8), (1 << 9), (1 << 10),
		(1 << 11), (1 << 12), (1 << 13), (1 << 14), (1 << 15), (1 << 16) };


	int nwm1[17] = { 0, (1 << 1)-1, (1 << 2)-1, (1 << 3)-1, (1<< 4)-1,
		(1 << 5)-1, (1 << 6)-1, (1 << 7)-1, (1 << 8)-1, (1 << 9)-1, (1 << 10)-1,
		(1 << 11)-1, (1 << 12)-1, (1 << 13)-1, (1 << 14)-1, (1 << 15)-1, (1 << 16)-1 };

	try {
		if (log_table[0] == NULL) log_table[0] = new uint16_t [nw[w]];
		if (ilog_table[0] == NULL) ilog_table[0] = new uint16_t [nw[w]*3];
	} catch (bad_alloc &e) {
		cerr << e.what() << endl;
		if (log_table[0] != NULL) delete log_table[0];
		if (ilog_table[0] != NULL) delete ilog_table[0];
		log_table[0] = NULL;
		ilog_table[0] = NULL;
		return -1;
	}

	for (j  = 0; j < nw[w]; j++) {
		log_table[0][j] = (uint16_t) nwm1[w];
		ilog_table[0][j] =  (uint16_t) 0;
	}

	b = 1;
	for (j = 0; j < nwm1[w]; j++) {
		if (log_table[0][b] != (uint16_t) nwm1[w]) {
			fprintf(stderr, "Galois_create_log_tables Error: j=%d, b=%d, B->J[b]=%d, J->B[j]=%d (0%o)\n",
					j, b, log_table[0][b], ilog_table[0][j], (b << 1) ^ prim_poly[w]);
			exit(1);
		}
		log_table[0][b] = (uint16_t) j;
		ilog_table[0][j] = (uint16_t) b;
		b <<= 1;
		if (b & nw[w]) b = (b ^ prim_poly[w]) & nwm1[w];
	}

	for (j = 0; j < nwm1[w]; j++) {
		ilog_table[0][j+nwm1[w]] = ilog_table[0][j];
		ilog_table[0][j+nwm1[w]*2] = ilog_table[0][j];
	}
	ilog_table[0] += nwm1[w];
	return 0;
}

int galois_create_log_tables(int w, uint32_t **log_table, uint32_t **ilog_table)
{
	int i, j, b;

	int prim_poly[33] =
	{ 0,
		/*  1 */     1,
		/*  2 */    07,
		/*  3 */    013,
		/*  4 */    023,
		/*  5 */    045,
		/*  6 */    0103,
		/*  7 */    0211,
		/*  8 */    0435,
		/*  9 */    01021,
		/* 10 */    02011,
		/* 11 */    04005,
		/* 12 */    010123,
		/* 13 */    020033,
		/* 14 */    042103,
		/* 15 */    0100003,
		/* 16 */    0210013,
		/* 17 */    0400011,
		/* 18 */    01000201,
		/* 19 */    02000047,
		/* 20 */    04000011,
		/* 21 */    010000005,
		/* 22 */    020000003,
		/* 23 */    040000041,
		/* 24 */    0100000207,
		/* 25 */    0200000011,
		/* 26 */    0400000107,
		/* 27 */    01000000047,
		/* 28 */    02000000011,
		/* 29 */    04000000005,
		/* 30 */    010040000007,
		/* 31 */    020000000011,
		/* 32 */    00020000007 };  /* Really 40020000007, but we're omitting the high order bit */

	int nw[33] = { 0, (1 << 1), (1 << 2), (1 << 3), (1 << 4),
		(1 << 5), (1 << 6), (1 << 7), (1 << 8), (1 << 9), (1 << 10),
		(1 << 11), (1 << 12), (1 << 13), (1 << 14), (1 << 15), (1 << 16),
		(1 << 17), (1 << 18), (1 << 19), (1 << 20), (1 << 21), (1 << 22),
		(1 << 23), (1 << 24), (1 << 25), (1 << 26), (1 << 27), (1 << 28),
		(1 << 29), (1 << 30), (1 << 31), -1 };


	int nwm1[33] = { 0, (1 << 1)-1, (1 << 2)-1, (1 << 3)-1, (1<< 4)-1,
		(1 << 5)-1, (1 << 6)-1, (1 << 7)-1, (1 << 8)-1, (1 << 9)-1, (1 << 10)-1,
		(1 << 11)-1, (1 << 12)-1, (1 << 13)-1, (1 << 14)-1, (1 << 15)-1, (1 << 16)-1,
		(1 << 17)-1, (1 << 18)-1, (1 << 19)-1, (1 << 20)-1, (1 << 21)-1, (1 << 22)-1,
		(1 << 23)-1, (1 << 24)-1, (1 << 25)-1, (1 << 26)-1, (1 << 27)-1, (1 << 28)-1,
		(1 << 29)-1, (1 << 30)-1, 0x7fffffff, 0xffffffff };

	try {
		if (log_table[0] == NULL) log_table[0] = new uint32_t [nw[w]];
		if (ilog_table[0] == NULL) ilog_table[0] = new uint32_t [nw[w]*3];
	} catch (bad_alloc &e) {
		cerr << e.what() << endl;
		if (log_table[0] != NULL) delete log_table[0];
		if (ilog_table[0] != NULL) delete ilog_table[0];
		log_table[0] = NULL;
		ilog_table[0] = NULL;
		return -1;
	}

	for (j  = 0; j < nw[w]; j++) {
		log_table[0][j] = (uint32_t) nwm1[w];
		ilog_table[0][j] =  (uint32_t) 0;
	}

	b = 1;
	for (j = 0; j < nwm1[w]; j++) {
		if (log_table[0][b] != (uint32_t) nwm1[w]) {
			fprintf(stderr, "Galois_create_log_tables Error: j=%d, b=%d, B->J[b]=%d, J->B[j]=%d (0%o)\n",
					j, b, log_table[0][b], ilog_table[0][j], (b << 1) ^ prim_poly[w]);
			exit(1);
		}
		log_table[0][b] = (uint32_t) j;
		ilog_table[0][j] = (uint32_t) b;
		b <<= 1;
		if (b & nw[w]) b = (b ^ prim_poly[w]) & nwm1[w];
	}

	for (j = 0; j < nwm1[w]; j++) {
		ilog_table[0][j+nwm1[w]] = ilog_table[0][j];
		ilog_table[0][j+nwm1[w]*2] = ilog_table[0][j];
	}
	ilog_table[0] += nwm1[w];
	return 0;
}

int galois_create_split_w08_tables(uint32_t **galois_split_w8)
{
	Galois *gc;
	int p1, p2, i, j, p1elt, p2elt, index, ishift, jshift;
	uint32_t *table;

	try {
		for (i = 0; i < 7; i++) {
			galois_split_w8[i] = new uint32_t [( 1 << 16 )];
		}
	} catch (bad_alloc &e) {
		for (i = 0; i < 7; i++) {
			if (galois_split_w8[i] != NULL) delete galois_split_w8[i];
			galois_split_w8[i] = NULL;
		}
		return -1;
	}

	gc = new Galois_Shift(32);
	for (i = 0; i < 4; i += 3) {
		ishift = i * 8;
		for (j = ((i == 0) ? 0 : 1) ; j < 4; j++) {
			jshift = j * 8;
			table = galois_split_w8[i+j];
			index = 0;
			for (p1 = 0; p1 < 256; p1++) {
				p1elt = (p1 << ishift);
				for (p2 = 0; p2 < 256; p2++) {
					p2elt = (p2 << jshift);
					table[index] = gc->galois_single_multiply(p1elt, p2elt);
					index++;
				}
			}
		}
	}
	delete gc;
	return 0;

}

void galois_region_xor(unsigned char *r1, 
		unsigned char *r2, 
		unsigned char *r3, int nbytes)
{
	/* nbytes must be a multiply of uint64_t (8 bytes) */
	uint64_t *l1;
	uint64_t *l2;
	uint64_t *l3;
	uint64_t *ltop;
	unsigned char *ctop;

	ctop = r1 + nbytes;
	ltop = (uint64_t *) ctop;
	l1 = (uint64_t *) r1;
	l2 = (uint64_t *) r2;
	l3 = (uint64_t *) r3;

	while (l1 < ltop) {
		*l3 = ((*l1)  ^ (*l2));
		l1++;
		l2++;
		l3++;
	}
	return;
}

void galois_invert_binary_matrix(int *mat, int *inv, int rows)
{
	int cols, i, j, k;
	int tmp;

	cols = rows;

	for (i = 0; i < rows; i++) inv[i] = (1 << i);

	/* First -- convert into upper triangular */

	for (i = 0; i < cols; i++) {

		/* Swap rows if we ave a zero i,i element.  If we can't swap, then the 
		   matrix was not invertible */

		if ((mat[i] & (1 << i)) == 0) {
			for (j = i+1; j < rows && (mat[j] & (1 << i)) == 0; j++) ;
			if (j == rows) {
				fprintf(stderr, "galois_invert_matrix: Matrix not invertible!!\n");
				exit(1);
			}
			tmp = mat[i]; mat[i] = mat[j]; mat[j] = tmp;
			tmp = inv[i]; inv[i] = inv[j]; inv[j] = tmp;
		}

		/* Now for each j>i, add A_ji*Ai to Aj */
		for (j = i+1; j != rows; j++) {
			if ((mat[j] & (1 << i)) != 0) {
				mat[j] ^= mat[i];
				inv[j] ^= inv[i];
			}
		}
	}

	/* Now the matrix is upper triangular.  Start at the top and multiply down */

	for (i = rows-1; i >= 0; i--) {
		for (j = 0; j < i; j++) {
			if (mat[j] & (1 << i)) {
				/*        mat[j] ^= mat[i]; */
				inv[j] ^= inv[i];
			}
		}
	}
}
void reed_sol_galois_w08_region_multby_2(unsigned char *region, int nbytes)
{
  uint32_t *l1;
  uint32_t *ltop;
  unsigned char *ctop;
  uint32_t tmp, tmp2;
  int prim08 = -1;
  int mask08_1 = -1;
  int mask08_2 = -1;
  Galois *gc;
  gc = galois_init(8, -1);


  if (prim08 == -1) {
    tmp = gc->galois_single_multiply((1<<7), 2);
    prim08 = 0;
    while (tmp != 0) {
      prim08 |= tmp;
      tmp = (tmp << 8);
    }
    tmp = (1 << 8) - 2;
    mask08_1 = 0;
    while (tmp != 0) {
      mask08_1 |= tmp;
      tmp = (tmp << 8);
    }
    tmp = (1 << 7);
    mask08_2 = 0;
    while (tmp != 0) {
      mask08_2 |= tmp;
      tmp = (tmp << 8);
    }
  }

  ctop = region + nbytes;
  ltop = (uint32_t *) ctop;
  l1 = (uint32_t *) region;

  while (l1 < ltop) {
    tmp = ((*l1) << 1) & mask08_1;
    tmp2 = (*l1) & mask08_2;
    tmp2 = ((tmp2 << 1) - (tmp2 >> 7));
    *l1 = (tmp ^ (tmp2 & prim08));
    l1++;
  }
}
void reed_sol_galois_w16_region_multby_2(unsigned char *region, int nbytes)
{
  uint32_t *l1;
  uint32_t *ltop;
  unsigned char *ctop;
  uint32_t tmp, tmp2;
  int prim16 = -1;
  int mask16_1 = -1;
  int mask16_2 = -1;
  Galois *gc;

  gc = galois_init(16, -1);

  if (prim16 == -1) {
    tmp = gc->galois_single_multiply((1 << 15), 2);
    prim16 = 0;
    while (tmp != 0) {
      prim16 |= tmp;
      tmp = (tmp << 16);
    }
    tmp = (1 << 16) - 2;
    mask16_1 = 0;
    while (tmp != 0) {
      mask16_1 |= tmp;
      tmp = (tmp << 16);
    }
    tmp = (1 << 15);
    mask16_2 = 0;
    while (tmp != 0) {
      mask16_2 |= tmp;
      tmp = (tmp << 16);
    }
  }

  ctop = region + nbytes;
  ltop = (uint32_t *) ctop;
  l1 = (uint32_t *) region;

  while (l1 < ltop) {
    tmp = ((*l1) << 1) & mask16_1;
    tmp2 = (*l1) & mask16_2;
    tmp2 = ((tmp2 << 1) - (tmp2 >> 15));
    *l1 = (tmp ^ (tmp2 & prim16));
    l1++;
  }

}
void reed_sol_galois_w32_region_multby_2(unsigned char *region, int nbytes)
{
  uint32_t *l1;
  uint32_t *ltop;
  unsigned char *ctop;
  int prim32 = -1;
  Galois *gc;
  gc = galois_init(32, -1);

  if (prim32 == -1) prim32 = gc->galois_single_multiply((1 << 31), 2);

  ctop = region + nbytes;
  ltop = (uint32_t *) ctop;
  l1 = (uint32_t *) region;

  while (l1 < ltop) {
    *l1 = ((*l1) << 1) ^ ((*l1 & 0x80000000) ? prim32 : 0);
    l1++;
  }
}

